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ABSTRACT 

In order to establish whether the unstable r-modes in a rotating neutron star provide a 
detectable source of gravitational waves, we need to understand the details of the many 
dissipative processes that tend to counteract the instability. It has been established 
that the bulk viscosity due to exotic particles, like hyperons, may be particularly 
important in this respect. However, the effects of hyperon superfluidity have so far 
not been fully accounted for. While the associated suppression of the reaction rates 
that give rise to the bulk viscosity has been estimated, superfluid aspects of the fluid 
dynamics have not been considered. In this paper we determine the r-mode instability 
window for a neutron star with a £~ hyperon core, using the appropriate multifluid 
formalism including, for the first time, the effect of the "superfluid" bulk viscosity 
coefficients. We demonstrate that, even though the extra terms may increase the bulk 
viscosity damping somewhat, their presence does not affect the qualitative features of 
the r-mode instability window. 



1 INTRODUCTION 

Comprising one and a half solar masses inside a radius of roughly ten kilometers, neutron stars provide an arena where 
many extremes of physics meet. A detailed model of neutron star dynamics must account for strong magnetic fields, various 
superfluid/superconducting components, the interaction between the crust nuclei and the fluid, as well as exotic states of 
matter that may be present in the neutron star core. Needless to say, it is a formidable task to construct such a model. 
Especially since it requires an understanding of physics well beyond the laboratory. While the equation of state for matter 
approaching the nuclear saturation density no w 0.16 fm~ 3 (corresponding to 2.48 x 10 14 g/cm 3 ) is quite well understood, it 
seems unlikely that laboratory experiments will ever be able to probe the densities expected in the deep core of a neutron 
star (above several times no). 

Despite decades of research into the supranuclear equation of state, considerable uncertainties remain. Furthermore, 
neutron star (NS) observations have only recently begun to reach the level of precision necessary to constrain the theoretical 
models in a severe way. Very recent results by Ozel, Baym & Oliver (2010) suggest that current measurements of NS masses 
and radii can be used to rule out several nuclear equations of state and support the notion that the core of a NS should contain 
exotic particles, such as pions, kaons, hyperons or even a deconfined quark condensate. Establishing observable signatures of 
the presence of such exotic states of matter is a priority for modelling in this area. 

During the last decade, the notion that gravitational waves (GWs) may drive the so-called r-modes of a rotating neutron 
star unstable, and that this may lead to the star spinning down on a timescale of weeks to months, has been discussed in 
a number of papers, see Andersson & Kokkotas (2001) and Andersson (2003) for reviews. The r-mode instability initially 
attracted attention because it provided a mechanism that could spin a newly born neutron star down dramatically, releasing 
GWs at a level that might be detectable in the process (Owen et al. 1998). For the purpose of GW detection, rapidly rotating 
accreting neutron stars in low-mass X-ray binaries (LMXBs) have also attracted attention. In these systems, the r-modes could 
provide a mechanism for torque balance (Bildsten 1998; Andersson et al. 1999) and, in some cases, lead to persistent GW 
emission (Andersson et al. 2002; Wagoner 2004; Nayyar & Owen 2006). This scenario currently holds interesting prospects 
for detection, albeit with considerable technical difficulties (Watts et al. 2008). 

Not surprisingly, the details of the supranuclear equation of state (EOS) are key to understanding the r-mode instability. 
In a real neutron star many (viscous) mechanisms compete with the GW driving of the r-mode. If the star contains exotic 
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Figure 1. The r-mode instability window for stars with M = IAMq and R = 13km. On the left we have a star where the dominant 
damping mechanism at low temperatures is shear due to electrons in an Ekman layer at the base of the crust, on the right a star where 
the main damping mechanism above 10 9 K is bulk viscosity due to hyperons. The rotation rate is expressed in terms of the parameter 
ui = u)/^/GM/R3. The instability curve in the left panel always has a negative slope and the system undergoes a limit cycle as described 
in the text. This is in contrast with the instability curve on the right which has a positive slope in the 10 9 K region, which could halt 
the thermal run-away and lead to persistent emission of gravitational waves. 

particles, such as hyperons, additional dissipation channels may become relevant. At first sight, it would appear that additional 
damping should reduce the chances of the r-mode GWs being detectable as it would reduce the region of parameter space 
where the instability is active. However, this is not necessarily the case. As an illustration of this, let us consider the instability 
in the temperature range 10 7 - 10 10 K. First assume that the main damping mechanism is due to a viscous Ekman layer at 
the core-crust interface (for a discussion see Glampedakis & Andersson (2006a,b)) (the left panel of figure 1). Standard bulk 
viscosity due to modified URCA reactions is only relevant above 10 10 K and is not included in this example (we shall discuss 
this in detail in section 6). In this case a NS in an LMXB, which will heat up to a core temperature of a few times 10 s K, 
would spin up due to accretion and enter the instability window before the star reaches the break-up limit. The shear from 
the unstable r-mode then heats the star and the emission of GWs spins the star down until it returns to the stable region. At 
this point the star will cool down and the cycle can begin again (Levin 1999). Unfortunately, with the current estimates for 
the nonlinear saturation amplitude (Arras et al. 2003; Brink et al. 2004) the duty cycle for this scenario is very low, meaning 
that the star would emit brief bursts of gravitational radiation and we would observe most systems in quiescence. Let us 
contrast this model with a system where we have added the effect of hyperon bulk viscosity (the right panel of figure 1). In 
this case, extra damping leads to a positive slope of the curve in the 10 9 K region and there are now three possible scenarios. 
Depending on the mode amplitude and on the exact details of the damping, the star could either i) execute the cycle we 
have already described, or ii) the heating may be sufficient for the system to evolve horizontally all they way to the positively 
sloped part of the curve before GW emission has time to spin the star down (Bondarescu et al. 2009). Finally, it may be the 
case that, iii) the heating due to accretion is such that the system becomes unstable in the region with positive slope. In the 
last two scenarios the system will not be able to evolve away (significantly) from the instability curve, and should become a 
persistent source of GWs (Andersson et al. 2002). This hypothesis was examined by Nayyar & Owen (2006), who found that, 
in the case of a NS with a hyperon core, persistent emission is possible over a wide range of parameters for the bulk viscosity 
and the EOS. 

LMXBs contain old NSs which will have cooled well below the temperature at which the hyperons (and other components 
of the star, such as neutrons and protons) are likely to become superfiuid. Superfluidity adds dimensions to the problem as, 
not only does it reduce the reaction rates for hyperon creation processes, it also increases the dynamical degrees of freedom of 
the system. In general, we have to work with multifluid hydrodynamics, with each particle species (potentially) leading to an 
independent flow. This leads to the appearance of new families of oscillation modes (Epstein 1988; Lindblom & Mendell 1994; 
Andersson & Comer 2001) and also has profound consequences for viscous dissipation. Not only are there new dissipation 
mechanisms, such as mutual friction between the components (Alpar, Langer & Sauls 1984; Mendell 1991a,b; Andersson et 
al. 2006), but bulk- and shear viscosity can no longer be described by single coefficients. In general, there are many additional 
viscosity coefficients. In the context of neutron stars, this was first pointed out by Andersson & Comer (2006). The particular 
case of a hyperon core was first considered by Gusakov & Kantor (2008) and has recently been analysed in detail by Haskell 
et al. ((in preparation)). In the simplest case, that of a core comprising neutrons, protons, electrons and E~ hyperons, one 
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Figure 2. We plot the rotation rate at which hyperons and dcconfined quarks first appear for neutron stars of different baryon masses 
{Mg). The data was obtained using the EOS G240 (Glcndcnning 1996) (left panel) and case III of Glcndcnning (1985) (right panel). 
For both equations of state stars with masses just above the threshold for the appearance of hyperons would have their exotic core spun 
out (or significantly reduced) as they approach the Kcplcrian rotation frequency. The horizontal line indicates the baryon mass that 
corresponds to a gravitational mass of M = IAMq in a non-rotating model. 

can show that the problem is very similar to that for superfluid Helium (Andersson & Comer 2008) and can be described by 
three bulk viscosity and one shear viscosity coefficient. One would, of course, also need to account for mutual friction between 
the various superfluid components. It has, however, been shown in several calculations that in a NS composed of neutrons, 
protons and electrons, mutual friction is unlikely to have a significant effect on the r-mode instability (Lindblom & Mendell 
2000; Lee & Yoshida 2003; Haskell et al. 2009). As the nature of mutual friction involving hyperons is largely unknown, we 
shall assume that it can also be neglected (the veracity of this assertion obviously needs to be checked by detailed work in the 
future). We will focus on the bulk viscosity, which is expected to give the main contribution to the r-mode damping (Jones 
2001; Lindblom & Owen 2002; Haensel et al. 2002; Nayyar & Owen 2006). Although the relevance of the new "superfluid" 
bulk viscosity coefficients is well established for superfluid Helium, their effect has mostly been neglected in the study of NS 
oscillations. The notable exception is the work of Kantor & Gusakov (2009) who studied the damping of sound waves in a 
dense superfluid hyperon core and showed that the additional bulk viscosity terms can play a significant role. 

It is clearly important to understand the role of the extra damping coefficients and refine our theoretical understanding 
of the bulk viscosity, as we have seen that the nature of the damping mechanisms can have profound consequences for the 
r-mode instability and the associated GW emission. In fact, direct GW detection from these systems should allow us to 
discern whether the star is emitting persistently or is executing a limit cycle. This would provide valuable information on 
the physics of the NS interior. The purpose of this paper is to study the effect of superfluid hyperon bulk viscosity on the 
r-mode instability window. The formalism for studying r-modes in multifluid neutron stars has been developed by Haskell et 
al. (2009) and we shall extend it to include hyperon bulk viscosity, thus considering for the first time the global dynamics of 
a multifluid NS with a hyperon core. 

2 SIZE OF THE HYPERON CORE 

Let us begin by considering the extent of the hyperon-rich core. It is well known that the central density of a star decreases 
as the rotation rate increases. This means that the threshold density for the presence of hyperons moves closer to the centre 
of the star and it could, in theory, be possible to "spin out" such an exotic core completely. Examples of this effect are shown 
in Figure 2. The figure shows how, for two relativistic mean field EOS, the size of the hyperon core (in the equatorial plane) 
depends on the star's rotation. Let us focus on stars with a constant baryon mass (horizontal lines in the figure) corresponding 
to a gravitational mass of IAMq in the non-rotating limit. Then, in the case of the G240 EOS (Glendenning 1996) (the left 
panel), there are no A hyperons present when the break-up limit is reached and the S~ core is very small. In the other 
example, for the EOS given by case 3 in Glendenning (1985), both the A and the E~ hyperons are completely spun out before 
the break-up limit is achieved. 

It is obviously relevant to quantify this effect and consider its impact on the r-mode instability window. Before doing 
this it is, however, important to discuss how reliable Newtonian estimates of the size of an exotic core are. This is relevant 
since r-mode viscous damping timescales are mostly calculated in Newtonian theory so far (although Nayyar & Owen (2006) 
calculate the background stellar model in general relativity and Pons et al. (2005) calculate shear viscosity damping for 
relativistic r-modes). A useful answer is provided in Figure 3. The results demonstrate that the relative size of the core is 
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Figure 3. Comparing results for neutron star models calculated from the EOS G 2 40 in both Newtonian gravity (dashed lines) and full 
General Relativity (solid lines). The first two panels illustrate the well-known fact that the two theories lead to rather different results 
for neutron star masses and radii given the same central density. From these panels it is clear why it is not particularly meaningful to 
use a realistic equation of state in a Newtonian study. The third panel is relevant for studies of effects due to the presence of a hypcron 
core. It is shown how the ratio of the hyperon core (here R c corresponds to the radius at which the density is such that the E~ hypcron 
first appears) to the radius of the star varies with the central density. Interestingly, these results indicate that, despite the stellar models 
being very different, the relative size of the core is almost the same in Newtonian theory and General Relativity. 



well approximated by the Newtonian model, even though the actual size of the star is off by a large amount compared to the 
relativistic result. 

This suggests that, even though it is generally not meaningful to use realistic EOS in Newtonian models (since the radius 
of a star with a given mass/dentral density differs so much from the relativistic model), one can study the relative size of the 
hyperon core also in Newtonian theory. To do this, we construct a sequence of stars spinning at different rotation rates but 
with the same total mass. For simplicity, we now restrict ourselves to rotating n = 1 polytropes (a useful model since both the 
background and the r-mode solution can be studied analytically). Following the analysis in Haskell et al. (2009) we consider 
a variable a which labels the (rotationally) deformed equipotential surfaces, and which is defined by the relation 

r = a[l + e(a,0)] (1) 

where e is a function which represents the rotational deformation of the equilibrium structure from the spherical background 
model. To second order in the slow-rotation approximation, the deformation can be cast in the form 

e = Di(o) + D 2 (a)P 2 (cos9) (2) 

where P 2 is the I — 2 Legendre polynomial. For an n = 1 polytrope we have (Chandrasekhar 1933) 

Di = — —4 u) and D 2 = =— ui (3) 

TV 2 RM r 9 RM r V ' 

where M r is the mass contained within a radius r, and R and Mo are the radius and mass of the non-rotating star, respectively. 
We have also defined 
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where we have used the dimensionless variable y = air/R. In terms of the variable a the equations of hydrostatic equilibrium 
take the simple form 

1 dP{a) _ d<f> R (a) 



p(a) da 



(5) 



with $h = $ — |fi 2 a 2 sin 2 8. Since the variable a is associated with equipotential surfaces of it follows that the density 
profile for an n — 1 polytrope has the same functional form as in the spherical case. That is, we have 

sin y tvM 

P = Pc— Wlth ^=4R3" ( 6 ) 



and we find that the mass of the rotating star is given by 

M — Mo 



, , 2 .2 

71-2 13 



(7) 



Using this relation, we can impose that M remain constant for all rotation rates and thus determine the central density as a 
function of a). This then determines the value of the coordinate a of the transition density where hyperons first appear. Figure 
4 shows an example of the extent of the exotic core, for different rotation rates, for this analytic polytropic model. 



Superfluid hyperon bulk viscosity and the r-mode instability of rotating neutron stars 5 



lr 

0.9 - 

0.8 ; 

0.7 - 



■ 0.5 - 



0.4 - 
0.3 : 



0.2 - 
0.1 - 
I — 



M=1.2Ms,R=14km 

M=1.4Ms,R=14km 

M=1.4Ms,R=10km 



Q. I ( GM/R 3 ) 



1/2 



Figure 4. Extent of the exotic core, for two models with M = IAMq, for R = 14 km and R = 10 km and a model with M = 1.2Mq 
and R = 14 km (the equation of state is taken to be an n = 1 polytrope). The graph extends to the Keplerian breakup frequency. For 
the stellar model with M = 1.4Mq and R = 14 km the extent of exotic core is reduced at the breakup frequency, while for the model 
with M = 1.2Mq the hyperon core is essentially spun out before the breakup rotation rate is reached. 



3 ESTIMATING THE HYPERON BULK VISCOSITY 

The aim of the present work is to estimate the effect that hyperon bulk viscosity has on the unstable r-modes. In the core 
of a mature neutron star, where we may expect a sizable hyperon population, several components are likely to be superfluid 
at temperatures below ~ 10 9 K. This complicates the picture considerably as it becomes necessary to use a multifluid 
description of the system. This, in turn, gives rise to, potentially, quite a large number of dissipation coefficients (Andersson 
& Comer 2006; Haskell et al. (in preparation)) All calculations to date have considered single fluid systems, including the 
effects of superfluidity only in the calculation of the reaction rates. This is somewhat restricted as it means that the only 
dissipation coefficients that have been considered are the "standard" ones for bulk and shear viscosity. A notable exception 
is the calculation of the "superfluid" bulk viscosity coefficients for a neutron star with a hyperon core by Gusakov & Kantor 
(2008), and the application of these results to the damping of sound waves (Kantor & Gusakov 2009). These studies show that 
the additional damping coefficients can play a significant role. In the following we shall estimate the superfluid bulk viscosity 
coefficients in the case of a core comprising neutrons, protons and E~ hyperons, a system that turns out to closely resemble 
superfluid Helium (Andersson & Comer 2008). 



3.1 Multifluid equations of motion 

Let us briefly outline the multifluid equations of motion for a system formed of neutrons (n) and protons (p) locked to the E~ 
hyperons. For simplicity, we neglect the presence of the neutral A hyperons (which, if they are superfluid, would add another 
degree of freedom) and only consider the bulk viscosity due to the process 

n + n ^ p + E" (8) 

This is, of course, not the only contribution to the hyperon bulk viscosity. However, reactions involving A hyperons have been 
considered by Lindblom & Owen (2002) and Gusakov & Kantor (2008). Their results demonstrate that the inclusion of A does 
not impact strongly on the qualitative features of the bulk viscosity. Our main reason for ignoring the presence of A hyperons 
is that we then have to consider only two fluids, the superfluid neutrons and a charge-neutral conglomerate of protons and 
E~ hyperons. In this case we can use the analytic results of Haensel et al. (2002) for the bulk viscosity coefficients. 

We shall assume that protons and E~ hyperons are locked together by the Coulomb interaction, which proceeds on a 
much faster timescale than the dynamical timescale we are interested in (i.e. the r-mode frequency). We therefore consider 
one single charge-neutral conglomerate. A full derivation and in detail discussion of the relevant equations of motion can be 
found in Haskell et al. ((in preparation)) and Andersson & Comer (2006). We will also ignore leptonic reactions, such as 
modified and direct URCA, as they proceed on a much longer timescale compared to the reaction in (8) (Haensel et al. 2002). 
Finally, we will (more or less) neglect the presence of electrons. They could easily be included by assuming that they are also 
locked to the protons (as in the usual model for the outer core of a neutron star (Mendell 1991a, b)), and that their inertia can 
be neglected (Andersson et al. 2010). In the presence of E~ the number density of electrons is also depleted (due to overall 
charge neutrality) , which means that they become less relevant in the deep core anyway. 
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Before locking the hyperons to the protons we have three distinct components. Following Andersson & Comer (2006) we 
can write the momentum for each of component as 

Tii = ga (m n n n « n - 2[a np w np + a nE w^ E ]) (9) 

< = 9H (m p n p ^ + 2[(a np + a pE K P - a P VE]) (10) 
= 9ij (m^n^v J s + 2[(a nE + a pS )w J nS - a pE ™ np ]) (11) 



In these expressions, n x (x = n, p, E) is the number density of each species, v x is the velocity of each component and 
!fx y = «x — Vy represents relative flows. The coefficients a xy describe the so-called "entrainment" effect, which leads to the 
momentum of a given species not being aligned with the individual velocity. The equations of motion then take the form of 
coupled Euler equations (omitting for the moment the effects of gravity); 



(12) 



where /i x is the chemical potential of each species, ff represents the sum of all external forces acting on the component x and 
D* J represents the dissipative part of the stress tensor, which includes shear and bulk viscosity (Andersson & Comer 2006) . 
In the following we focus on hyperon bulk viscosity and thus only explicitly include the corresponding terms in the equations 
of motion. These equations are complemented by continuity equations for each component, in the form: 

<9 t n x + V t ( n x u x ) = T x (13) 



where T x is the particle creation rate per unit volume for component x. In general, T x depends on all reactions involving the 
x particle species. This will render the complete expressions quite complicated. This is another reason why we focus on the 
single reaction (8). Finally we assume overall baryon number conservation, i.e. 

E r * = ° ( 14 ) 



3.2 The linearised problem 

Since we are interested in the r-modes we will follow Haskell et al. (2009) and consider perturbations of a rotating background 
in which all the fluids flow together. Assuming that the fluids are in chemical equilibrium we then have F x = in the 
background. In writing the perturbation equations we shall assume, as discussed previously, that the £~ hyperons are also 
locked to the charged component by the Coulomb interaction. In this case one has (indicating Eulerian perturbations with S 
and Lagrangian perturbations with A) Sw' nS = <5ui np and the problem is reduced to that of a two-fluid flow. The equations of 
motion can be written in terms of two independent velocities which, following Haskell et al. (2009) we take to be 5w' np and 
the "total" velocity Sv t , defined as 

x 

where p x = m x n x is the mass density of each component and p = p x . We also introduce the total pressure p such that 

ViP^^rixVi^x (16) 

x 

In terms of these variables the perturbed Euler equations (including the dissipative terms due to bulk viscosity), to linear 
order and in a frame rotating with the star, can be cast in the form: 

dt(pSvi) + V t 5 P + 2pe ljk WSv k - ^V,p = -V 3 ^ Df^j (17) 

d t [(1 - e-e 2 )5w^] + \vS + 2e ijk W5w k np = -V J f ^ - ^ - ^ ] (18) 
2 I p n 2ps 2p p / 

and the continuity equations take the form 

d t Sn b = -Vj (n b Sv j ) - ^Vij* (19) 

d t Sx s = — (l + x s — ) Vif ~ Sv l \7 t x s + ^ (20) 
rib \ m n J n b 
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where we have defined 

f = n b x s (l - yc)Sw l np , (21) 
Am = m s — m n , (22) 

e = ^±^1, (23) 
P2/c(l - Vc) 

Am a np - a nS 

e 2 = , (24) 

m n py c 

/3 = 2/i n /m n - ns/ms - p- P /m p . (25) 

We have also introduced the baryon number density rib, the fraction x x = n x /i%b (with x c = (tie + n p )/rib) and the mass 
fraction y c = (pa + p P )/p. Note that /3 does not vanish in the background. Rather, it is of order 0(Am/m n ), as chemical 
equilibrium with respect to the reaction in (8) implies 2/u n — jUe — fi p = 0. If we consider bulk viscosity to be the only dissipative 
process at work, we can write the dissipative contributions to the stress tensor as 

£>« = X)^ = ~9ii [CViW + C nB V ( /] (26) 



" ' " - -^[^VjW+^Vd'] (27) 



D n D E D p 

p n 2p E 2p p 

This allows us to cast the Euler equations in the form 

d t { P Svi) + V,Sp + 2pe ijk n j 5v k - ^ V t p = V, [(V,^ + C nE V ; j ; )] (28) 

d t [(1 - e - £2)8w^} + X - + 2e l3fc n^^ p = V, [c s V,Sv l + C^a] (29) 

As we can see, we now have four bulk viscosity coefficients, of which only three are independent (as expected from the analysis 
in Haskell et al. ((in preparation))). That is, we need to determine two bulk viscosity coeffients that are not present in the 
single fluid problem (only the ("-term is present in the Navier-Stokes equations). These extra coefficients can be calculated 
from the equation of state and the reaction rate Fe- 

Finally, we impose charge neutrality. As we are neglecting the electrons in the core, we have for the background 

x s w x p w y (30) 
Meanwhile, for the perturbations charge neutrality leads to the condition 

At 

A^ p = A^ E = (31) 

where A represents a Lagrangian perturbation. It is defined by 

A = S + C ( (32) 

where £j is the Lie derivative with respect to the Lagrangian displacement £ l associated with the co-moving motion, such 
that dtC — Av' (note that as we are in a two-fluid system it would also be possible to define a Lagrangian displacement 
associated with the counter-moving motion). 

Solving the full set of equations (19)-(29), including the dissipative terms, is still a prohibitive task. We thus follow 
the common strategy of assuming that the dissipative terms only introduce a small deviation from the solution of the non- 
disipative problem, and use this solution to estimate the bulk viscosity damping timescale. The non-dissipative equations 
follow if we set Fe = C = C" E = C nE — C E = 0- We also note that, since T x = 0, we can rewrite (20) as 

d t 5y s = X - (\ + V,- {py c [(l - y.)Swi p ]} - 6v*V jys . (33) 

These are, in fact, almost exactly the equations that were considered by Haskell et al. (2009). The only difference is the 
entrainment dependence in equation (29), where one has the extra term £2 due to the difference in mass between neutrons 
and S~ hyperons. The similarity of the final equations obviously means that it is straightforward to adapt the method from 
Haskell et al. (2009) to the present problem. 



4 SINGLE-FLUID BULK VISCOSITY REVISITED 

Given our aim, we take as starting point the study of Haensel et al. (2002). Their results are particularly useful because they 
provide relatively simple parameterised expressions for the S~ hyperon bulk viscosity coefficient also in the case of superfluid 
constituents. Of course, these explicit expressions come at a price. They are based on the bare-particle assumption, e.g. do not 
consider the effective masses (of neutrons, protons and hyperons). As emphasised by Haensel et al, this approximation is likely 
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to be rather severe and one may find that "dressed particle" effects affect the results significantly. In order to allow for this 
possibility, Haensel et al opted to leave a free parameter (x) in their formulae. In our study, we will use the freedom associated 
with this parameter to discuss the plausible range for the hyperon bulk viscosity. This "phenomenological" approach to the 
problem is reasonable given the many uncertainties associated with the supranuclear equation of state. Finally, let us remark 
that the calculation of the bulk viscosity coefficient is made assuming that the fluids are locked together. As discussed in 
section 5.1, this is equivalent to neglecting the dependance of the reaction rate on the divergence of the relative velocities. 
Although this assumption is not physically justified, we make it for the sake of simplicity and in order to make progress on 
the r-mode problem. 

To calculate the bulk viscosity coefficient due to the reaction n + n ^ p + £~ we follow the procedure of (Sawyer 1989). 
In the single-fluid case the reaction rate Fs depends only on the lag in instantaneous chemical potentials 

A/3 = (2A/u n - Afj, p - Afis) (34) 

Then 

T E = — A A/3 (35) 

where the coefficient A can be obtained from Haensel et al. (2002). If we assume that the fluid is in equilibrium with respect 
to all other reactions, and impose that the perturbations maintain charge neutrality, i.e. 

An E = An p , (36) 

we can write the lag in chemical potentials as a function of two parameters, e.g. the total baryon number and the hyperon 
fraction. Thus, A/3 = A/3(n b ,x s ), from which we obtain 

A/3 = BAx-s + C— (37) 
n b 

with 

B=^- and C=^-n h (38) 
ctes drib 

By inserting equation (37) into the continuity equations (19)-(20) and making use of the fact that Ax c = 2Axs, which 
follows from equation (31), we have 

d t An b = -n b Vi (W) (39) 

(40) 
(41) 



dtAxv = -A— (bAxv + C— ) 
n b \ n b J 



If we assume a harmonic time dependance with frequency ui for the perturbations, these relations allow us to compute the 
hyperon fraction. Finally, as the dissipation is given by the part of the pressure due to deviations from chemical equilibrium, 
we have (denoting with p eq the pressure in chemical equilibrium): (note that the variation is taken at constant baryon number 
density) : 

P = Pe q + ^-Ax s (42) 

OXs 

which, inserted into the Euler equations (28) and (29) leads to the result of Haensel et al. (2002); 

C 2 n 2 b 1 



^ |A|B 2 l + a 2 

where 

um b 



(43) 
(44) 



|A|S 

Here |A| oc 1/T (in non-superfluid matter) is related to the relaxation time-scale of the relevant reactions which drive matter 
towards equilibrium. We can then use the estimate of Haensel et al. (2002) : 



a ~ 6 09 ( m " ] 2 mp ms 104 ( nb ) ( - — 1 ( i-vw-m i 

\m*J m* r?4 T 9 2 x V 1 far 3 / \ n s J \ B J ( ' 

where T 9 = T/10 9 K and m = w/10 4 s~\ In the following, we shall take as canonical values = 0.7m x and x — 0.1, 
as in Haensel et al. (2002). We will also focus on a stellar model represented by an n — 1 polytrope. The reason for this 
simplification is that we can make direct use of the analytical r-mode results of Haskell et al. (2009). In future studies it 
would, of course, be desirable to calculate the coefficients in (43) consistently from a realistic equation of state and carry out 
the mode-analysis for models constructed using the same equation of state. However, in order for this level of modelling to be 
meaningful one also has to account for general relativity, c.f. the discussion in Section 2. In the case of the r-mode instability 



1 far 3 \ 1/3 /l00MeV\ 
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this has been done, but not for truly realistic equations of state (Ruoff & Kokkotas 2001, 2002; Lockitch et al. 2001, 2003; 
Yoshida et al. 2005; Gaertig & Kokkotas 2008). Further work is clearly needed. 

Before we proceed, let us note that the bulk viscosity is weak whenever the parameter a is either very small or very large. 
That is, when the relaxation timescale associated with the reaction (8) is significantly shorter or longer than the oscillation 
timescale. This is natural since bulk viscosity is essentially a resonance effect associated with the dynamics and the reactions. 
As a result, at constant baryon number density rib and oscillation frequency uj, there are two asymptotic regimes; 




at high temperatures, and 

cr - e.3 x W W (iSEv) i- ■ 

at low temperatures. Note that we could equivalently refer to the asymptotic expressions as low- and high-frequency approx- 
imations (at fixed temperature). We expect these asymptotic expressions to be good at low temperatures (below 10 9 K) and 
high temperatures (above 10 10 K), but they are not valid in a temperature range where 

Jow ^ >high , rr ^r„,„9 -1/2 -1/6/ n b X 1 / 3 / w X 1 / 2 , . 

In this region the relaxation timescale is comparable to the oscillation timescale, and one would not expect the approximations 
above to be useful. 



5 SUPERFLUIDITY 

At temperatures below ~ 10 10 K hyperons (and also neutrons and protons) are expected to form Cooper pairs and become 
superfluid. As we have already discussed, this can have a significant impact on the reaction rates and the bulk viscosity. The 
critical temperature below which each species is superfluid depends strongly on the density. In figure 5 we show the critical 
temperature for neutrons (using model "e" from Andersson et al. (2005)), protons (model "h" from Andersson et al. (2005)) 
and S~ hyperons, where we have assumed, as in Nayyar & Owen (2006) that the pairing gap is the same as that for A 
hyperons computed by Balberg & Barnea (1998). This is essentially the gap deduced from BCS theory, so it does not account 
for medium effects etcetera. This is worth keeping in mind since such effects are known to reduce the pairing gaps for neutrons 
and protons significantly. From the figure we see that there will be superfluid hyperons present below roughly 5 x 10 9 K. We 
also learn that the transition will not take place until much lower temperatures are reached in regions above 4-5 times nuclear 
density (above n b ~ 0.6 fm~ 3 ). In fact, the figure indicates that the hyperons in the deep core may not be superfluid in most 
astrophysical systems, e.g. the accreting neutron stars in LMXBs which are expected to have core temperatures above 10 8 K. 
Note that the pairing gaps we have used represent rough approximations to the detailed many-body results and different 
models lead to gap results that are generally within factors of a few of one another (for a review see Andersson et al. (2005)). 
This will have an effect on the suppression of the reaction rates and the extent of the different regions described below, but 
it is unlikely to affect the qualitative picture that emerges from our calculations. 

When one or more of the species involved in the reactions that lead to the macroscopic bulk viscosity are superfluid, the 
reaction rates will be suppressed. This can be accounted for by a multiplicative factor |A| — > 1Z\\\. The form of this suppression 
factor depends strongly on which species involved in reaction (8) are superfluid and thus on the temperature and the region 
of the star one is considering. In the region where only the hyperons are superfluid (region IV of figure 5) we can use the 
analytic reduction factor of Haensel et al. (2002), valid for singlet state pairing (note that we use the notation of Haensel et 
al. (2002) and the symbols are not to be confused with those used in other sections of the paper): 

5/4 , ,1/2 

7l s = ± exp(0.5068 - \/0.5068 2 + y 2 ) (49) 



y = (l.456 - + ) with r - (50) 



where 

0.157 1.764 \ . , T 

—— + With T = — 

V"T T J T c 

while a = 1 + 0.3118 y 2 and 6= 1 + 2.556 y 2 and T c is the superfluid critical temperature for hyperons, which can be related 
to the superfluid gap by the relation y — A(T)/kbT, where kb is Boltzmann's constant. Note that the above expression would 
be valid also in the case where only protons are superfluid and give us 1Z P , as long as one uses the corresponding critical 
temperature. For the neutrons in the core, on the other hand, there will be triplet pairing, so to describe the reduction due 
to neutron superfluidity we shall consider the approximate factor of Haensel et al. (2002) 

Tin = (0.6192 + x/0.3808 2 + 1.156V) x exp(0.7756 - \/0.7756 2 + y' 2 ) 

+0.18766 y 2 x exp(1.7755 - \/l.7755 2 + 4y 2 ) (51) 
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Radius (Km) 

Figure 5. Critical temperature for superfluidity of neutrons, protons and E — hyperons for a star with M=1.4 Mq and R=10 km. For 
the hyperons we use the model of Balberg & Barnea (1998) while for the neutrons we use model "e" and for the protons model "h" of 
Andersson et al. (2005). In regions I and II the neutrons are superfluid, but not the hyperons, and the system should exhibit multifluid 
behaviour. Region III, in which all components are superfluid, will also exhibit multifluid behaviour. In region IV only the hyperons are 
superfluid and all components should flow together. 



To describe the reduction rate lZ np needed in regions I and II of figure 5 we would need to account for the reduction due to 
the superfluidity of both neutrons and protons. Haensel et al. (2002) do not provide an analytic expression for the reduction 
rate in these regions, but Nayyar & Owen (2006) found that the product of two single particle reduction factors (i.e. the 
factors calculated for the case when only one species is superfluid) is a good approximation for the reduction rate in the case 
when both species are superfluid. In particular, they found that one could take lZ p s ss TZ^TL^. Assuming that this result is 
generally valid, we use 

K„ p w TtrJlp (52) 

Note that, in the r-mode problem, the prescription for the regions where only neutrons and protons are superfluid (regions I 
and II of figure 5) is not crucial, as these regions do not dominate the contribution to the hyperon bulk viscosity damping. 
Region I is relatively unimportant because the r-mode eigenfunctions increase with radius, while the contribution from region 
II is negligible as the hyperon number density falls steeply below rib ~ 0.32 fm~ 3 for our model EOS. However, we need a 
prescription for the region where all the constituents are superfluid (region III of figure 5). In this region the reduction factor 
of Haensel et al. (2002) can only be evaluated numerically, so we shall use the prescription 

^npE « K n H p 1l s (53) 

which reproduces the qualitative features of the result of Haensel et al. (2002). In particular, the suppression becomes stronger 
when all the species are superfluid. 

As a means of comparison we shall examine the effect of our prescription on the r-mode instability. To do this we shall 
calculate the standard single-fluid r-mode solution, as described by Haskell et al. (2009), and calculate the critical frequency 
at which hyperon bulk viscosity stops the mode from growing unstable (cf. section 6 for more details). In figure 6 we plot the 
critical r-mode frequency curve and compare the result obtained with our prescription for the superfluid suppression to that 
obtained with the approximate reduction rate used by Lindblom & Owen (2002) 

— = e ' (54) 

where A is the hyperon pairing gap. The comparison shows that the approximation (54) is reasonably accurate, although it 
is also clear that the more detailed reduction rate in (49) leads to a slightly smaller reduction factor. 

To complete our description we also need to specify in which region to use the one-fluid bulk viscosity timescale of equation 
(71) and in which to use the two-fluid timescale of equation (74). Technically, as the hyperons are locked to the protons, the 
system will exhibit two-fluid behaviour whenever the neutrons are superfluid (regions I, II and III of figure 5). However, as 
we have mentioned above, we will only consider regions I and III, and neglect region II as it is essentially hyperon-free and 
does not contribute to the bulk viscosity. The whole star could of course be included and the picture be made consistent by 
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Figure 6. The r-mode instability window in the single fluid case including the effect of hyperon bulk viscosity. We compare the superfluid 
reduction rate used in this work with the approximation of Lindblom & Owen (2002) (LO in the figure). It is clear that the approximation 
of Lindblom & Owen (2002) is reasonably accurate at lower temperatures when all species are superfluid. In the higher temperature 
region we use the exact reduction rate of Haensel et al. (2002) for the case where only hyperons are superfluid, which leads to a smaller 
reduction factor and thus a slightly smaller instability window. 

using the full results for the superfluid reduction rates of Haensel et al. (2002), but such a calculation is beyond the scope of 
this qualitative study. 

5.1 Superfluid bulk viscosity 

To obtain a complete description of the bulk viscosity we need to have an estimate not only of the "usual" bulk viscosity 
coefficient, but also of the two extra "superfluid" ones. We will now use the standard bulk viscosity coefficient of Haensel et al. 
(2002), discussed in the previous section, to estimate the superfluid bulk viscosity coefficients for the reaction n + n ^ p + S~. 
This strategy obviously has some flaws, as the reaction rates and coefficients of Haensel et al. (2002) are calculated in the one 
fluid approximation, in which the reaction rate Fe depends only on the instantaneous lag between the chemical potentials. 
For our multifluid system the reaction rate should depend, in general, also on the divergence of the relative velocity. That is, 
it should take the form (Andersson & Comer 2006) 

r E = -\ w A/3 - T TO ViA< E (55) 

In the following we shall assume that X m = A (i.e. the same as in the one fluid case) and that t w = 0. In practice, we 
are assuming that the variation of the reaction rate with respect to the single fluid result is small. This assumption is not 
necessarily justified but as there is, to the best of our knowledge, presently no calculation of the coefficients A^ and t w it 
is the only option if we are to make progress. For future applications it would be highly desirable to have a calculation of 
the reaction rates for the full superfluid system and estimates of the impact of the parameter t w on the r-mode damping 
timescale. Given the reaction rate Fs, we follow the strategy of the previous section and expand the pressure and difference 
in chemical potentials around a background solution at equilibrium. The continuity equations now have an extra term due to 
the difference in velocity between the neutrons and the charged component 

d t An b = —n b \7 i Av l ~^^'Vj\n b xs(l — y c )Awi v \ (56) 

m n L 

d t Ax s = (BAx s + C^) + \n b x s (l - y c )AwU (57) 

rib \ n b J n b 1 

Note that we have expressed the continuity equations in terms of Lagrangian perturbations, but as we are working in the 
rotating frame and have chosen a co-moving background, we have Aui„p = Sw^p and Av l = Sv 1 . Following the analysis of the 
previous section we can expand the parts of the pressure and /3 which depend on the deviations from chemical equilibrium; 

Sp = 5p eq - (VM - C nS Vii' (58) 

and 

= (H)"[CV^ + C nE V^] (59) 
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where we recall that f — n b xs{l — y c )Sw l np , and £ is the same coefficient as in equation (43). Using equation (59) in the Euler 
equations (29) we can infer the following relation 

C nE C 

V- = -V (60) 

c s C 

which, as expected from the Onsager symmetry principle (Haskell et al. (in preparation)), shows that in fact only three of the 
coefficients are independent. Moreover, we can rewrite equation (59) as 



/3 = /3 eg - 

which allows us to identify the coefficients: 



C ViSv ViJ +C ViJ 

m n n b 



rib 



B 

C 



( Am\ Am] - nS / Am\ Am 

\ m n / m n J \ m n / m n nf, 



c 2 



l + x E - 



Am 



(61) 



(62) 



(63) 



If we identify £ with the coefficient obtained by Haensel et al. (2002), i.e. the one given in (43), the above relation allows 
us to determine the other coefficients from the equation of state, and evaluate the bulk viscosity damping timescale. Basically, 
this is equivalent to assuming that the reaction rates are the same as in the single fluid case. The coefficients needed in 
equation (29) are then: 

'dp 



JnE _ 1 

^ ~ 2 



(64) 



- E= l 
s 2 



m n n b 



(65) 



Finally, we can compare our results to those obtained from the relativistic formulation of Gusakov & Kantor (2008). This is 
clearly a rough comparison, as we have neglected both relativistic effects and the effect of A hyperons but, if we for simplicity 
neglect entrainment and terms of the order of Am/m n , we can identify the dissipative stress tensor Dij in equation (26) 
with the relativistic analogue in equation (63) of Gusakov & Kantor (2008). Comparing the two expressions we find that the 
various coefficients agree reasonably well. 



6 THE R-MODE INSTABILITY WINDOW 

The stability of an r-mode can be determined by estimating the gravitational-wave driving and viscous damping timescales. 
The result is usually illustrated in terms of the critical frequency at which the driving and damping timescales are equal. In 
our case, the instability curve is obtained by solving for the roots of 

— + — + — + — = (66) 

T gw THb T sv TEk 

where THb is the hyperon bulk viscosity damping timescale, t 3V is the shear viscosity damping timescale, teh is the damping 
timescale due to an Ekman layer at the base of the crust and r gw is the gravitational-wave growth timescale, which for an 
n = 1 polytrope and the I = m = 2 r-mode, is given by (Andersson & Kokkotas 2001): 

M \ _1 / R V 4 ( p 



^ = - 47 U4MeJ UkmJ llmsj S (67) 
(the sign implies that the mode is growing). The nature of the shear viscosity damping will depend on which region of the 
star we are considering, and one should also note that in a multifluid system there will, in general, be more shear viscosity 
terms than in the single-fluid problem. The exact nature of the scattering processes that give rise to shear viscosity depends 
on which species are present and whether they are superfiuid. Shear viscosity in a hyperon core has not, to the best of our 
knowledge, been studied in the literature. Hence, we cannot at this stage consider hyperon scattering processes. We could, 
however, quantify the effect that the presence of hyperons has on the standard scattering processes. In particular, we expect the 
electrons to be severely depleted in the presence of E~, thus weakening the shear from electron-electron and proton-electron 
scattering. This should reduce the shear viscosity in the part of the star where neutrons, protons and E~ are all superfiuid. In 
regions where the neutron are normal, e.g. the deep core, we also know that neutron-neutron scattering will dominate and we 
can use the estimates of Andersson et al. (2005). These details may not be that relevant, however, since the various scattering 
processes lead to weaker damping than the shear associated with the Ekman layer at the crust-core boundary (Bildsten & 
Ushomirsky 2000; Lindblom et al. 2000; Levin & Ushomirsky 2001). Since there will be no electron depletion in the outer 
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Figure 7. We compare the critical curves for shear viscosity due to n-n scattering and for shear due to an Ekman layer, as described in 
the text. The damping from the Ekman layer is always stronger in the region of interest. For reference we have also included the critical 
timescale for hyperon bulk viscosity (with suppression of the reaction rates included, following the prescription described in the text). 



core, we can use the standard result for the boundary layer damping. Thus we take 

2 / n \ 1/2 



We arrive at this estimate by taking the simple constant density estimate of Andersson & Kokkotas (2001) for a M — IAMq 
neutron star, corrected for a "slippage" factor S c =0.05, as defined by Glampedakis & Andersson (2006b). It has been shown 
by Glampedakis & Andersson (2006a) that one should expect the constant density estimate to only differ by factors of a 
few from the result for a stratified model. Hence, it should be a reasonable approximation. Having said that, it is absolutely 
clear that the boundary layer issue needs more detailed scrutiny. Based on our current understanding, the physics in the 
crust-core transition region dictates the r-mode damping at the temperatures that are relevant for mature neutron stars. Yet, 
our understanding of the effect that superfluidity and superconductivity may have on the boundary layer is far from complete 
(Kinney & Mendell 2003; Sidery 2008). These issues are of central importance to neutron star dynamics, and more detailed 
studies should be encouraged. 

The hyperon bulk viscosity damping timescale is given by 

1 1 (dE s 

nTb ~ ~2E [~dt 

To leading order in rotation, the energy of the mode takes the form (Haskell et al. 2009) 

E k = ~Jp [Sv 2 + (1 - g)i, c (l - y c )Sw 2 np ] dV (70) 

where we recall that y c — y p + and the definition (23). In the one fluid region (which also includes the deep core), where 
the neutrons flow together with the protons and hyperons, we can write 

dE 



(69) 



at 



({v l 8vydv (7i) 



where the bulk viscosity coefficient £ is taken to be that of Haensel et al. (2002) in the hyperon core, while in the outer layers 
of the neutron star we take the value given by Sawyer (1989), appropriate for modified URCA reactions in neutron, proton 
and electron (npe) matter 



^ 6x1 ° 2 Hh^J 2 («) 6 ^ 2 (72) 

where uj r is the r-mode frequency in the rotating frame. This leads to the r-mode damping timescale [where we have corrected 
an error in the numerical prefactor of the result from Andersson & Kokkotas (2001)] 



^ 8 - 6xl0l \i^JUkW feJ IotJ s (73) 

We note that although we have added the npe bulk viscosity contribution for completeness it is, in fact, irrelevant in the 
temperature range of interest for mature neutron stars. This is evident from the results in figure 8. The result is natural given 



M W R Y 1 ( P \ 2 ( T 
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Temperature (K) Temperature (K) 

Figure 8. In the left panel we show the effect on r-mode instability window of including npe bulk viscosity in the outer layer of the 
star. As is clear from the figure the effect is completely irrelevant below Fa 10 10 K, above which the contribution due to modified URCA 
reactions becomes dominant. This can easily be understood from the right hand side panel where we compare the damping timescale for 
hyperon bulk viscosity and for npe bulk viscosity. It is obvious from the figure that the npe bulk viscosity damping timescale is always 
much longer in the region of interest (below fa 10 10 K). 



the temperature scaling in (72) and the fact that the resonance associated with the modified URCA reactions is located at 
higher temperatures (above 10 10 K) compared to the hyperon resonance (at a few times 10 9 K) (note, however, that if a 
sizable proton fraction is present in the core, then direct URCA reactions can contribute to the bulk viscosity in the region 
where all baryons are superfluid and the hyperon reactions are reduced (Haensel et al. 2002)). 

Let us now consider the two fluid region, where one has a host of extra terms due to the separate fluxes of the neutrons 
and the charged particles. In this region we need to evaluate 



dE 

m 



J C (W) 2 + ic nE [(V^ 1 ) (Vj*) +c.c] + 

ffi {C E [(v<j<) (Vj*) + c.c] + C nE [(Vj l ) (V l 5<) + c.c] } dV 



+ . 



+m 



V / c 



(JyJ KAi E )* + c.c. 



+ 



c 1 



Am \2 



[(pAy^)(n b Ax^)* + c.c] 



c 



ii E 



i/.Af/E) ( ^ ] + r.r. 



dV (74) 



where c.c indicates the complex conjugate, 2j ! = py c (l — yc)&Wn P and f(e) = (1 — s)/(l — e — As). Note that, following 
Haskell et al. (2009) we have assumed that <5w" p vanishes identically outside the two fluid region and that Ap = Axe = 
at the surface of the star, which is equivalent to all the components having the same surface. This condition is "natural" 
as we do not expect the outer regions of the crust to be superfluid. We have taken e to be a constant, although this is not 
necessarily, realistic, and in the following we will assume that, as an approximation to the results of Gusakov et al. (2009a), 
a nc = a nS + C(Am/m n ) and thus /(e) = 1 + C(Am/m n ) 2 . This is clearly an approximation, that follows from assuming 
(as described in appendix A) that the coefficients of the relativistic entrainment matrix are constant. Use of the full density 
entrainment parameters from Gusakov et al. (2009a), without neglecting higher order terms in Am/m n , would produce a 
weak entrainment dependence of the damping timescale. 



6.1 Model EOS 

In order to assess the superfluid hyperon bulk viscosity effect on the r-modes, we will make use of the analytic solution by 
Haskell et al. (2009). This means that we assume that the overall density profile is that of an n — 1 polytrope, in which case 
the solution for the comoving degree of freedom in the superfluid problem is completely independent from the countermoving 
motion (which is, in fact, driven by the comoving part). In order to evaluate the integral in (74) we need to consider the speed 
of sound, which for an n — 1 polytrope of the form p — Kp 2 is 

c 2 = 2Kp (75) 
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We shall also assume, as an "approximation" to case III EOS of Glendenning (1985), that in the core the hyperon number 
density follows a relation of the form 

xt — qin b + q 2 (76) 

with qi — 0.114 fm 3 and q 2 = —0.036. This allows us to consider a simple model that is linear in the number density of 
baryons, with the key feature that there are no hyperons below n b « 0.319 fm 3 . A more realistic model would give a steep 
drop in the number density of hyperons close to the transition density, but as we are not using a realistic equation of state 
our linear approximation is sufficient. With the use of the thermodynamical identity 

dp 



2 UP 

dx-s dn b 



one can then obtain 
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c 2 s m n 
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( , 9x ^ 
xt. + n b - — 
V on b 


<tin 2 b 


m n 
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where we have used the relation p ■ 

8xt 



m n n. b (l + isAm/mn). We can now simplify the integral in (74) making use of: 



Axt = 



'-An b 



(77) 

(78) 
(79) 

(80) 



Furthermore, we can write the mass fraction j/e as 

rriT 



m n 



ie (1 — k) and p = m n n b (l + n) 



where we expand in terms of the parameter k = Am/m^XT- This is justified as k 
From (81) we can then conclude that 

niT 



Finally, charge neutrality (31) implies p p ■■ 



which leads to 



and 



Aj/ E 

Pt and from (16) 
'dp 



Axt. (1 - 2k) 



(81) 

0.03 for the densities we shall consider. 



(82) 



1 

m n 
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m n (1 — x 

m E x c 



8xt n b 



An b _ Ap 

n b PI dn b xt 
Using the values of the bulk viscosity coefficients from equation (63) we can now write: 

BE 2 f /Ap\ 



dt 



with 



Ceff 



1 + 



Am(! 



Jv 2 



6§2£n 6 ) 



C = 4 



1 + 



Am (1 — x c — 6an b ) 



(83) 



(84) 



(85) 



(86) 



(87) 



The integral in (86) is considerably simpler than the original expression in (74). In particular, it only involves calculating 
the integral of the quantity for the r-mode. As discussed above, and presented in detail in Haskell et al. (2009), the 
co-moving degrees of freedom decouple form the counter-moving ones to second order in rotation, so only a solution to the 
single fluid r-mode problem is needed. In particular, as we are considering Ann — 1 polytrope we shall use the explicit analytic 
solution in section 4.6 of Haskell et al. (2009). Note that the co-moving r-mode solution sources the counter-moving motion, 
thus giving rise to a relative flow and the additional bulk viscosity coefficients. 



7 RESULTS 

Before moving on to estimate the effect of the superfluid bulk viscosity, let us remark on the effect that rotation has on the 
instability window. In figure 9 we show an example of the r-mode instability window, in the one fluid limit, for a stellar model 
with M = 1.4M and R = 12.5 km, using the hyperon bulk viscosity coefficient of Haensel et al. (2002). We compare the 
critical frequency with rotational corrections included in the size of the core and the frequency obtained without including 
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Figure 9. The critical curve due to hyperon bulk viscosity in the single fluid case, for a star with M = XAMq and R = 12.5 km with a 
comparison between the non-rotating case and the case with the rotational corrections included. The straight horizontal line represents 
the Kcplcrian breakup velocity. 



them. It is immediately obvious that there is a significant difference for large rotation rates; in particular there is a range 
of temperature where the instability would be completely suppressed if we considered the non-rotating result, while this is 
not the case if we include the rotationally corrected size of the core. It is thus crucial to include the effects of rotation in the 
bulk viscosity calculation. This was already noted by Nayyar & Owen (2006). Note that the effect of rotation on the size of 
the core is technically an order 0(Q 2 ) effect compared to the leading order term (which would give the extent of the core 
in a non rotating star). This means that, to be consistent when calculating the bulk viscosity dissipation integral we would 
have to consider higher order terms also in the mode eigenfunctions, i.e. compute the whole integral to order 0(Q 4 ). This is, 
of course, prohibitive. We do, however, feel that it is appropriate to calculate the bulk viscosity damping for a sequence of 
stars with different rotation rates and that, when the core becomes very small (i.e. for large rotation rates) the effect of the 
reduced size will dominate (in other words we assume that the eigenfunctions of the mode are regular and that they will not 
contain terms that, at second order in rotation, become very large as the core becomes very small). It is important to bear in 
mind, however, that we are neglecting terms that could be of the same order as the change in size of the core and that could, 
especially for low rotation rates, lead to quantitative (but probably not qualitative) differences in the damping timescale. 

We can now consider the effect of introducing the extra "superfiuid" bulk viscosity coefficients in the two-fluid region. 
In figure 10 we show the instability window for two stellar models, both in the one-fluid approximation and including the 
two-fluid fluid bulk viscosity coefficients, while in figure 11 we show the effect of varying the stellar mass. Two effects are 
immediately obvious. First of all, the effect of the extra coefficients is to increase the bulk viscosity damping and reduce the 
instability window. Although the effect is not very pronounced in general, it is clear from the left panel of figure 10 that, by 
increasing the strength of the damping in the low frequency region of the curve, multifluid effects could possibly suppress the 
instability in certain systems. It is, however, also clear that the qualitative features of the instability window are, in general, 
not affected and there is still a rise in the critical curve in the T = 10 9 K region. This is due to the fact that the resonant 
nature of the bulk viscosity is robust. The result was expected, since Haskell et al. (2009) have shown that the r-mode is 
exclusively co-moving at the leading order in rotation. This means that the main conclusion of Andersson et al. (2002) and 
Nayyar & Owen (2006) is unaffected by the introduction of two-fluid effects. Notably, as the curve still has a positive slope in 
the 10 9 K region, the possibility still exists that for certain systems hyperon bulk viscosity could halt the thermal runaway of 
the r-mode and lead to a neutron star with a hyperon core becoming a persistent source of gravitational waves. Note, however, 
that if most of the hyperon core is not superfiuid then the system may cool rapidly [see eg. Page, Geppert & Weber (2006)]. In 
principle, this may halt a thermal runaway before the system reaches the positive sloping curve, leading once again to a limit 
cycle (Bondarescu et al. 2009). The key point is that observations may be able to distinguish systems evolving according to 
the distinct scenarios, thereby providing insight into the qualitative nature of the r-mode instability curve, and by inference 
possibly information about the state of matter in the deep neutron star core. 

Finally, before moving on, let us comment on the freedom associated with the parameter \- I n the left panel of figure 11 
we illustrate the effect of varying this parameter in the range 0.001 < % < 1- Clearly the effect on the height (frequency) of 
the peak of the resonance can be quite drastic, even though the qualitative nature of the curve is essentially unaffected, as it 
still exhibits a rise in the T = 10 9 K region. However for large values of \ the effect of hyperon bulk viscosity is weakened, 
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Figure 10. R-mode instability window, for a star with M = IAMq and varying radius. We show a comparison between the one-fluid case 
and the case with the two-fluid corrections included. The graphs extend to the breakup frequency and rotational corrections are included 
in both cases. The new "superfluid" coefficients increase the viscous damping but do not, in general, alter the qualitative features. 

making it much less likely that the thermal run-away of the system could be halted. We can gain an understanding of this 
from the right panel of figure 11, where we show the bulk viscosity coefficient £ for varying x a t nt — 1 fm~ 3 . It is clear that 
larger values of \ shift the peak to lower temperatures, where bulk viscosity is not the main damping mechanism, and make 
the effect weaker in the higher temperature region. Such a large range of values for \ is, however, unlikely and we choose to 
use the fiducial value of x — 0.1, as Gusakov & Kantor (2008) show that it reproduces their results and those of Lindblom & 
Owen (2002) to within factors of a few. It would obviously be good if future work were to better quantify the effects of the 
bare-particle assumption and narrow down the possible range for \- 



8 CONCLUDING REMARKS 

Our results represent the first investigation into the effect of including the extra superfluid bulk viscosity coefficients in a 
calculation of the r-mode instability window. We have shown that, even though the additional bulk viscosity coefficients do 
not alter the qualitative aspects of the instability window, there are regions of parameter space in which they could play 
a significant role, and may even suppress the instability entirely. In the light of these results we believe it is important to 
move beyond the qualitative analysis presented here. One should clearly account for the presence of A hyperons and use a 
more realistic equation of state to describe the star. Furthermore, if one is to construct a more realistic model one clearly 
needs to work in general relativity and include finite temperature effects (such as in the entrainment coefficients calculated by 
Gusakov et al (2009b)) and dissipative effects such as hyperon bulk viscosity. More detailed theoretical input is also needed 
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Figure 11. In the left panel we show the r-mode instability window for various values of the parameter \ in the hyperon bulk viscosity, 
for a M = XAMq and R = 13 km star. The graphs extend to the breakup frequency. We assume that the viscosity at low temperatures 
is due to an Ekman layer at the base of the crust and that the reaction rates are reduced by superfluidity as prescribed in section 5. 
We see that varying \ can have an significant effect on the height of the resonance but not on the qualitative nature of the curve. On 
the right we show the bulk viscosity coefficient £ as a function of temperature for 71^=1 fm — 3 and an oscillation frequency of 2300 Hz 
(chosen to directly compare with Nayyar & Owen (2006)). One can easily appreciate that the maximum of the coefficient is shifted by 
changing \, leading to the differences seen in the left panel, as explained in the main text. 



from the nuclear physics community, in order to calculate the superfluid reaction rates needed to evaluate the bulk viscosity 
coefficients. 

Developing the relevant tools will allow us to make progress on a range of related problems, e.g. involving finite temperature 
superfluid in the outer neutron star core or exotic phases of deconfined quarks in the deep core. An improved understanding 
of these systems is crucial if we are to take advantage of the unique opportunity that gravitational-wave detection could offer 
for the study of matter under extreme conditions. 
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APPENDIX A: THE ENTRAINMENT MATRIX 

In this Appendix we describe how to translate the results of Gusakov et al. (2009a) into our formalism. First of all we shall 
write the momentum of each fluid in the form 7tf = p xy vf, where the entrainment matrix p xy encodes the fact that the 
momentum of each component is not necessarily parallel to it's velocity in a superfluid. For a system of n, p, A and E~ the 
momenta take the form: 

(Al) 
(A2) 
(A3) 
(A4) 

Were the entrainment matrix p xy is connected to the relativistic entrainment matrix V xy of Gusakov et al. (2009a) by 

xy 

F xy = (A5) 
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By comparing with the momenta in equations (9)-(ll) one the finds that: 



2i/nn 

c Y 



\ ( m n n n - 2a np - 2a nE - 2a nA ) (A6) 
m n m p m n m s m n mjy 

2 v pp _ J- / _ , np „ pS r. pA> 



c"Y pp = m n n n - 2a np - 2a p2j - 2a pyl (A8) 

wip V / 

c 2 yP s = 2a pE ^ c 2y P A = 2a pA c 2ySA = 2q sa 

m p ms ' m p mA ' m^niA 



2 V SS 



2,/- A A 

c y 



^2- (m E n E - 2a nE - 2a pE - 2a EA ) (A10) 
^ (m A n A - 2a EA - 2a pA - 2a nA ) (All) 



where the a xy are the entrainment parameters that enter the equations of motion in section 3.1. We shall restrict ourselves 
to the case of no A hyperons and assume that the components of the relativistic entrainment matrix are constant in the core. 
In particular, as an approximation to the results in figure 1 of Gusakov et al. (2009a), we shall assume that V~ nE = Y np , from 
which we can obtain that: 

a np = a nE + C>(Am/m n ) (A12) 
which leads to the result used to calculate the dissipation integral in (74): 

/(e) = 1 + C(Am/m n ) 2 (A13) 
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